AN ANALYSIS OF THE QUASI-NONLOCAL QUASICONTINUUM 
APPROXIMATION OF THE EMBEDDED ATOM MODEL 



XINGJIE HELEN LI AND MITCHELL LUSKIN 

Abstract. The quasi-nonlocal quasicontinuum method (QNL) is a consistent hybrid coupling 
method for atomistic and continuum models. Embedded atom models are empirical many-body 
potentials that are widely used for FCC metals such as copper and aluminum. In this paper, we 
consider the QNL method for EAM potentials, and we give a stability and error analysis for a chain 
with next-nearest neighbor interactions. We identify conditions for the pair potential, electron 
density function, and embedding function so that the lattice stability of the atomistic and the 
EAM-QNL models are asymptotically equal. 



1. Introduction 

Hybrid atomistic-to-continuum methods couple atomistic regions surrounding defects with con- 
tinuum regions to achieve the accuracy of the atomistic model and the efficiency of the continuum 
model. Quasicontinuum hybrid methods utilize the Cauchy-Born rule for the energy density in 
the continuum region [21] . The original quasicontinuum energy [21] (denoted QCE) has interfacial 
forces (called "ghost forces") for a uniform strain 14,24]. Thus, uniform strain is not an equilibrium 
solution for the QCE energy (even though uniform strain is an equilibrium for purely atomistic and 
for purely coarse-grained continuum models). 

More accurate atomistic-to-continuum coupling methods have been proposed to remedy the 
QCE model. The ghost force correction method (GFC) achieves an increased accuracy by adding 
a correction to the ghost forces as a dead load during a quasistatic process [H0[lOl[18l[23]. The 
GFC method can be viewed as a stationary iterative method [HElDIllIIE] to solve the force-based 
quasicontinuum aproximation (QCF) using QCE as a preconditioner. More accurate coupling can 
be achieved by using a more accurate preconditioner or by using GMRES acceleration to solve the 
QCF equilibrium equations [SIIIIIEIIIIIj but the non-conservative and indefinite QCF equilibrium 
equations make the iterative solution and the determination of lattice stability more challenging [10 1 . 

An alternative approach is to develop a quasicontinuum energy that is more accurate than QCE. 
We will call a QC energy consistent if it does not have ghost forces for a uniformly strained lattice. 
The quasi-nonlocal energy (QNL) was the first consistent quasicontinuum energy [25] . For a one 
dimensional chain, the original QNL method is restricted to next-nearest neighbor interactions [25J. 
The QNL method for pair interaction potentials was extended to finite range interactions in [15] 
and to two dimensional finite range problems in [23] , 
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In this paper, we formulate a one-dimensional QNL energy for the embedded atom model (EAM) 
following |25j. The embedded atom model [12j[TJ:,20j is an empirical many-body potential that is 
widely used to model FCC metals such as copper and aluminum. We then give an analysis of 
the stability and error for the EAM-QNL approximation in the next-nearest neighbor case for a 
periodic chain. 

We identify conditions for the pair potential, electron density function, and embedding function 
so that the lattice stability of the atomistic and the EAM-QNL models are asymptotically equal. 
We also show in Remark 14.41 that the atomistic and EAM-QNL models can be less stable than the 
local quasicontinuum model (EAM-QCL), which is the EAM-QNL model with no atomistic region, 
if the above conditions on the pair potential, electron density function, and embedding function 
are not satisfied. 

Many theoretical analyses of QC models have been given based on pair-potential interactions 
[I T l6 ll9lfTT1ll5pi6pi9j . In this paper, we give an analysis of the stability and accuracy of a linearization 
of the quasi-nonlocal method for the EAM potential in one dimension with next-nearest neighbor 
interactions. A nonlinear a priori and a posteriori error analysis for the QNL model with next- 
nearest-neighbor pair potential interaction in one dimension was given in [22]. We think that a 
similar nonlinear analysis using the inverse function theorem can be done for this model, but we 
restrict our presentation in this paper to the linear analysis for simplicity. 

In Section [2j we present the notation used in this paper. We define the displacement space U and 
the deformation space 3^f- We then introduce the norms we will use to estimate the modeling error 
and the displacement gradient error. In Section EJ we introduce the QNL model with next-nearest 
neighbor interaction for the EAM potential. 

In Section HJ we give sharp stability estimates for both the fully atomistic model and the EAM- 
QNL model for a uniformly strained chain. Sharp stability estimates are necessary to determine 
whether quasicontinuum methods (or other coupling methods) are accurate near instabilities such 
as defect formation or crack propagation (7|9] . Similar stability estimates for the fully atomistic and 
fully local quasi-continuum (QCL) models can also be obtained by discrete Fourier analysis |13| . 

In section^ we study the convergence rate of the EAM-QNL model. We compare the equilibrium 
solution of the EAM-QNL model with that of the fully atomistic model, and we use the negative 
norm estimation method [EJ[l5] to obtain an optimal rate of convergence of the strain error. The 
error estimate depends only on the smoothness of the strain in the continuum region and holds 
near lattice instabilities, thus demonstrating that the QNL method for the EAM potential can give 
a small error if defects are captured in the atomistic region. 

2. Notation 

In this section, we present the notation used in this paper. We define the scaled reference lattice 

eZ := {e£ : £ G Z}, 

where e > scales the reference atomic spacing and Z is the set of integers. We then deform the 
reference lattice eZ uniformly into the lattice 

FeZ := {Fet : £ G Z} 

where F > is the macroscopic deformation gradient, and we define the corresponding deformation 
YF by 

{yF)t '■= f° r - oo < £ < oo. 
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For simplicity, we consider the space U of 2iV-periodic zero mean displacements u = (ui)i e z from 
y F given by 

hi := < u : ue + 2N = ue for i G Z, and = >, 

and we thus admit deformations y from the space 



y*F := {y : y = YF + u for some u G ZY}. 

We set e = 1/N throughout so that the reference length of the periodic domain is fixed. 
We define the discrete differentiation operator, Du, on periodic displacements by 

,„ x U£- U£-l 

(Du)£ := , — oo < i < co. 

We note that (Du) e is also 2A r -periodic in t and satisfies the zero mean condition. We will denote 
(Du) e by Du£. We then define 

Due - Du £ _i 



and we define (D^u) £ and (£> (4) u) £ in a similar way. To make the formulas concise and more 

readable, we sometimes denote Du£ by u' e , D^U£ by u", etc., when there is no confusion in the 
expressions. 

For a displacement ueW and its discrete derivatives, we define the discrete t\ norms by 

(at \ V 2 / n \ 1 / 2 

e X] l U£ n ' ll u/ H^ := e ' etc - 

£=-N+l / V £=-JV+l / 

Finally, for smooth real- valued functions £ (y) defined for y G 3^f 5 we define the first and second 
derivatives (variations) by 

" d£ 



(<5£(y),w) := ^ Q—(y)we for all w G W 



(5 2 £(y)v,w) := ^ — — — (y)v e w m for all v, w G W. 

£,m=-AT+l ^ ym 

3. The Embedded Atom Model and Its QNL Approximation 
We first give a description of the next-nearest neighbor EAM Model. 

3.1. The Next-Nearest-Neighbor Embedded Atom Model. The total energy per period of 
the next-nearest neighbor EAM model is 

£&(y) := £ a (y) + JF(y) (3.1) 

for deformations y G y F where £ a (y) is the total atomistic energy and J-(y) is the total external 
potential energy. The total atomistic energy is the sum of the embedding energy, £ a (y), and the 
pair potential energy, £ a {y) : 

£ a (y):=£ a (y)+£ a (y). (3.2) 
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The embedding energy is 

N 

£%y):=e £ G(p}(y)) 
e=-N+i 

where G(p) is the embedding energy function, the total electron density pf(y) at atom £ is 

Pe(y) ■= p(Ve) + PiVe + 1/1-1) + PiVe+i) + Piv'e+i + 1/1+2)1 

and p{r/e) is the electron density contributed by an atom at distance r. The pair potential energy 
is 

- 1 

£ a {y)-=e Y, 2 + ^ + + + ^+1 + ^+2)] 

£=-Af+l 

where ecj)(r/e) is the pair potential interaction energy [12]. Our formulation allows general nonlinear 
external potential energies .F(y) defined for y G 3^f ; but we note that the total external potential 
energy for periodic dead loads f is given by 

N 

f(y) ■= - Y e h ye - 

e=-N+i 

The equilibrium solution y a of the EAM atomistic model (|3.ip then satisfies 

-(5S a (y a ), w) = -(6£ a (y a ),w) - (S£ a (y a ), w) = (£F(y a ), w) for all w G W. (3.3) 

Here the negative of the embedding force of (|3.3j) is given by 

iV 

(J£ a (y a ),w)=e J] ^(^(y^^-^^^K + p^^ + ^JK + ^.J 
e=-N+i 

+ p'M+iK+l + p'Pyf+i + Dy1 +2 )(w' e+l + ^ +2 ) 
the negative of the pair potential force of fj3.3j) is given by 

N 

2 

e=-N+i 

+ (j>'(Dyf +1 )w' i+1 + 0'(£>% a +1 + £>% a +2 )K + i + w' e+2 ) 
and the external force is given by 

N dF 

(5J r (y),w) = Y g—{y) w l forallwGW. 
e=-N+i ye 

3.2. The Quasi-Nonlocal EAM Approximation for Next-Nearest-Neighbor Interactions. 

Hybrid atomistic-to-continuum methods can give an accurate and efficient solution if the deforma- 
tion y 6 yp is "smooth" in most of the computational domain, but not in the remaining domain 
where defects occur [U[22]. The goal of QC methods is to decompose the reference lattice into 
an atomistic region with defects and a continuum region with long-range elastic effects. It applies 
an atomistic model to the atomistic region for accuracy and a continuum model to the continuum 
region for efficiency. 

In this paper, we will consider an atomistic region defined by the atoms with reference positions 
x e for£ = —K, . . . , K, and a continuum region for t G {-N+l, . . . ,-(K+3)}U{(K+3), . . . ,N}. To 
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eliminate the ghost force that energy-based quasicontinuum approximations can have 
we define the remaining atoms, ±(K + 1),±(K + 2), to be quasi-nonlocal atoms [6l[25]. For the 
pair potential energy, the quasi-nonlocal atoms ±(K + 1), ±(K + 2) interact without approximation 
with atoms in the atomistic region, but interact through the continuum Cauchy-Born approximation 
with all other atoms [25j. The interactions of the quasi-nonlocal atoms for the embedding energy 
is slightly more complex, as given in [25j and below. 

The atomistic energy associated with each atom is given by 

£?(y) := £t(y) + £?(y) = G (#(y)) + \ [</>(%) + M + vU) + M+i) + M+i + lM 

where £f(y) denotes the embedding energy at atom £ and (y) denotes the pair potential energy 
at atom £ (£|(y), £f nl (y), £f{y) and £f n \y) will be defined analogously below), and the continuum 
energy associated with each atom is given by 

£ c t {y) :=£i(y)+£tty) =\G{p'i{y))) + \G{p c i+1 {y)) 

+ \ Wt) + <f>(H) + M+i) + <K%in)] 

where the total continuum electron density at atom £ is 

pj(y) := 2p(y' e )+2p(2y' e ). 

To define the QNL energy for the quasi-nonlocal atoms, we define the QNL electron density at 
atom £ by 

pf\y) :=2p(y' i )+2p(y' i + y' t _ 1 ). 
We then define the QNL energy for the quasi-nonlocal atoms by 

£t +l {y)-. = £f +l {y) + £fUy) 

= \G(pf +1 {y)) +\G(? K+ 2{y)) 

+ \ [Hv'k+i) + 4>(y'K+2) + Hv'k+i + v'k) + <f>(2y'K+2)] 



and 



< 2 (y): = 4 n j 2 (y) + ^t 2 (y) 

~~ 2 

1 
2 



G(pS 2 (y))+^(^ +3 (y) 

+ o [^(y'K+2) + Hv'k+s) + Hv'k+2 + VK+x) + 0( 2 fe+3) 



We define the QNL energy in a symmetric way and so only give the formulas for < £ < N. 
The total energy per period of the QNL model is then given by 

N 

e £ £f(y)+Hy) 

= £ <^(y) + J-(y) = £ <^(y) + £« nl {y) + ^(y), 
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where 

£1{y) for 0<£<K + 1, 
£ qnl {y) for e = K + l,K + 2, 
£f{y) for K + 2 < t < N. 

The equilibrium solution y qnl of the EAM-QNL model (|3.4p then satisfies 

- (y^ ) , w) = - (££« n ' (y qnl ) , w) - (<5£ gni (y 9 ™' ) , w) = (6 F{y qnl ) , w) for all w € W , 

where the negative of the embedding force is given by 

r qnl^ 
£=0 

+P / (Dy q e r l[)w' e+1 + p'(Dy q e l[ + Dy q ^ 2 ){w' f+l + wJ +2 ) 
+ ^'(^(y^)) • [p'pyflx^+i + p'P^j, + Dy? + <) 
+ eC(^ +a (y^)) • [p\Dyf +2 ) W ' K+2 + 2p\2Dyf +2 ){w' K+2 ) 

+ eG"(p^ 2 ( y ^)) • + pWkU + ^0«+2 + v/k+i 

+ eG'(ff K+3 (y qnl )) • [p , ( J D 2 /f| 3 X + 3 + 2p / (25yfi3)« + 3) 

AT 



+ e £ {G'(p-Ky 9 " i ))-[p / ( J DyrK + 2p'(2I) y f)K; 

+G'(^ +1 (y^)) • [p'{Dyf^) W [ +1 + 2p' (2Dyf^){w[ +1 )\ } , 



and the negative of the pair potential force is given by 
(S£ qnl (y qnl ),w) = ... 

+ eJ2\ [KDvFM + 4>'{Dyf + Dyf\){ W \ + vt t _ x ) 

£=0 

+4'(Dyf* 1 )v/ e+1 + fiiDyj+x + Dy q l l 2 )(w' e+1 + w' l+2 ) 
</>' (Dy q K l +1 )w' K+1 + 4>' \Dy q K +1 + Dyp)(w' K+1 + u^-) 
+ £ [<fi ' {Dyf +2 )w' K+2 + 2^(2D»g. 2 )(t4 +2 y 

4>' {Dy q K +2 )w' K+2 + 0'p2/F+2 + D y q K+l)( w 'K+2 + WK+i) 



+ 



+ 



2 L 



+ 



A? 
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4. Stability Analysis of The Atomistic and EAM-QNL Models 

In this section, we will give a stability analysis for the atomistic model and the EAM-QNL model 
for the next-nearest neighbor case. We will use techniques similar to those presented in [9] for the 
atomistic and QNL method for pair potentials. 



4.1. The Atomistic Model. The uniform deformation is an equilibrium of the atomistic 
model (|3.2p . therefore, we say that the equilibrium y_p is stable in the atomistic model if and only 
if (5 2 £ a (y F ) is positive definite, that is, 

(5 2 £ a (y F )u, u) = (8 2 £ a (y F )u, u) + (8 2 £ a (y F )u, u) > for all u € U \ {0}. (4.1) 

Note that {5 2 £ a (y F )u, u) is given by formula (7) in [9]: 

{5 2 £ a (y F )u,u) = iF||£>u||| - ey 2 ' F ||£> (2) u|||, (4.2) 

where 

Ap := <t>" F + m F for <t>" F := 4>"{F) and <f>% F := 4>"{2F) (4.3) 

is the continuum elastic modulus for the pair interaction potential. Thus, we only need to focus on 
(5 2 £ a (y F )u,\i), that is, 



(5 2 £ a (y F )u,u) = e 



£=-N+l 



Yl \ G "f [p'f{ u 'i + u 'i+i) + P2f{u£-i +u' e + u' e+1 + u' e+2 )] 2 

+G' F [p" F (u'd 2 + + uU? + P'M + i? (44) 

+P2 F {u' l+1 +u' l+2 ) 2 } 



where 



p'f 



p'(F), p' F :=p"{F), p' 2F :=p(2F), p'> F := p"(2F), 
G' F := G>{P}{y F )) = G'(p c e (y F )) = G'(pf(y F )), 



G" F := G"(p a e (y F )) = G"(p c e (y F )) = G"(pf(y F )). 

We calculate the identities 

(u' e + u' e+1 ) 2 = 2 (u' e ) 2 + 2 {u' e+1 ) 2 - e 2 « +1 ) 2 , 

(u' e + u' e+1 + u' e+2 ) 2 = 3 {u' t f + 3 {u' e+1 ) 2 + 3 (u' £+2 ) 2 - 3e 2 « +1 ) 2 - 3e 2 « +2 ) 2 + e 4 
2 (u' £ + u' i+1 ) • (v! t _ x +u' e + u' i+1 + u' e+2 ) 

{u'^f + 3 (u e ) 2 + 3 {u e+1 ) 2 + (u e+2 ) 



(4.5) 



t (3) 



«) 2 + 2« +1 ) 2 + (u" +2 y 



u e+i) + \ u e+2 
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We can now calculate explicitly the first equality below and then use (|4.5p (with u' replaced by u") 
for the second equality to obtain 



{u' e + u' e+1 + u e+2 + u £+3 ) 2 = 4 ((u' e ) 2 + (u e+1 ) 2 + (u e+2 ) 2 + (u £+3 f) 



e iu 



i+2) 



e [u 



+ Vt+2) 



2 I n i // 2 I ii . // , // V 

4(K) 2 + K +1 ) 2 + K +2 ) 2 + K +3 ) 2 ) 



(6« +1 ) 2 + 8« +2 ) 2 + 6« 4 



+ 6 4 (4(4 3 + ) 2 ) 2 + 4( 



(3) n2 
£+3 ) 



e\u\%) 2 . 



We can then obtain from the above identities that 

(6 2 £ a (y F )u, u) = G" F ■ { [4 (p' F ) 2 + 16 (p 2F ) 2 + 16p^] ||Du||| 

-6 2 [(p' F ) 2 + 20 (p 2F ) 2 + l2p' F p' 2F ] ||L>^u||| 

8 (// 2F ) 2 + 2^] ||L>( 3 )u||| - e 6 (p' 2F ) 2 pWu|||} 
3 ' 2 V) pu|||- 2e 2 p 2 VP (2) u|||} 



+ 6^ 



|2 



{4G'; (p' F + 2p' 2F ) 2 + 2G' F {p" F + 4p' 2 V) } \\Du\ 
- e 2 [G" F [( P ' f ) 2 + 20 (p' 2F ) 2 + 12p' F p' 2F ] + G' F 2p'i F } \\D®u\\% 



+ e A G F 



8(p' 2F ) 2 + 2p' F p' 2F \ \\D®u\ 



'G' F (p' 2F ) 2 \\D^u\ 



We define the continuum elastic modulus for the embedding energy to be 



A F := AG" F (p' F + 2p' 2F f + 2G' F {p" F + 4p 2 ' F ) 



and 



Af := A F + A F , B F := - [$ F + G F ((p F ) 2 + 20(p' 2F ) 2 + \2p' F p' 2F ) + G' F {2p'i F )) , 
C F :=G F (8(p 2F ) 2 + 2p F p 2F ), and D F := -G" F {p' 2F ) 2 . 

Then (14. ip becomes 

(5 2 £ a (y F )u,u) =A F \\Du\\ 2 i + e 2 B F \\D^u\\ 2 i + e 4 C F ||Z)( 3 )u||| +e 6 D F \\D^u\\%. 

We will analyze the stability of (<5 2 £ a (y F )u, u) by using the Fourier representation |13j 

N 



k=-N+l V v * 7 



(4.6) 



(4.7) 



(4.8) 
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It then follows from the discrete orthogonality of the Fourier basis that 

- f 

(5 2 £ a (y F )u,u) = \ck\ 2 -<A F + B F 4sin 2 



-N+l 



kir 



+ C F 



4 sin 2 



kir 
2N 



2N 



+ D F 



(4.9) 



4 sin 2 



kir 



We then see from (|4.9p that the eigenvalues A& for k = 
respect to the ||L>u||^2 norm are given by 

Afc = Air(s fc ) for s k = 4 sin 2 



2N t 

N + 1, . . . , JV of (5 2 £ a {y F )u, u) with 
kir \ 



2N J 



where 

X F (s) := A F + Bps + C F s 2 + Dps 3 . 
From the pair interaction potential, electron density function, and embedding energy function 
given in Figure 2 in [12J, we assume that 

4>" F > 0, (f>2p < 0; p'p < 0, p' 2F < 0; p' F > 0, p'^p > 0; and G" F > 0. 

We then have from the assumption (I4.10p that 

C F > 0, Dp < 0, and 8\D F \ < C F . 

We can check that (jlTil]) implies that \D F s\ < 4\D F \ < C F /2, for < s < 4, so 

XUs) 



Bp + 2C F s + 3D F s 2 > B F + -ys for all < s < 4. 



(4.10) 
(4.11) 

(4.12) 



We conclude from (|4.12p that the condition Bp > or equivalently 



J 2F + G'f 



(p'p) 2 + 20 (p' 2F ) 2 + \2p'pp' 2F + G'p 2p'i F = -Bp < 0, 



(4.13) 



and the assumptions (14.10j) imply that A(s) is increasing for < s < 4. We thus have the sharp 
stability result 

(4.14) 



{S z £ a (y F )u,u) > A F (si)||L>u|| 2 2 > (A F +A F \ ||£>u||^ for all u£U. 
We summarize this result in the following theorem: 

Theorem 4.1. Suppose that the hypotheses (|4.10p and (|4.13p hold. Then the uniform deformation 
y F is stable for the atomistic model if and only if 



Af(si) = A F + B F 
= A F + A F 



4 sin 



2 (— ) 

\2NJ 



+ c F 



+ D F 



2F 



[p'p) 2 + 20 {p 2F ) 2 + l2p'pp' 2 p 



+ G 



'p 2p 2F ^ 



+ 4 2 sin 4 (^) G'p [r, (p' 2F ) 2 + 2p' F p' 2F ] - 4 3 sin 6 G" F (p' 2F ) 2 > 0. 

Remark 4.1. The role of the assumption (]4.13p is to guarantee that u' t = sin(e^7r) is the eigenfunction 
corresponding to the smallest eigenvalue of (5 2 £ a (y F )u, u) with respect to the norm ||Du||^2. In 
fact, we can see from the above Fourier analysis that u'^ = sin(e£-7r) is not the smallest eigenvalue 
of (5 2 £ a (y F )u, u) with respect to the norm ||Du||£2 for sufficiently large N if (I4.13P does not hold 
since then A'(0) < 0. 
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The assumption (|4. 13|) on the the pair interaction potential, electron density function, and em- 
bedding energy function cannot be expected to generally hold for physical embedded atom models 
since the nearest neighbor term G F (p' F ) 2 > dominates. We note, however, that generally G' F < 
for F > 1 [20], in which case G' F Ifl^p < 0; so ()4.13|) is more likely to hold for tensile strains F > 1. 



4.2. The EAM-QNL Model. Now we will analyze the stability of the EAM-QNL model for 
next-nearest neighbor interactions. The Fourier techniques used to analyze the stability of the 
atomistic model cannot be used for the EAM-QNL model because the Fourier modes are no longer 
eigenfunctions. Recall that the total atomistic interaction energy of the QNL model is £ qnl (y) := 

£ qnl (y) + £ qnl {y) = e J2e=- N +i where ^(y) is symmetric in £ e {-N + 1, . . . , TV} and is 

given by 



£T\y) ■= { 



£f{y) for 0<£<K + 1, 
£f +1 {y) for £ = K + 1, 
£f +2 {y) for £ = K + 2, 
{ £f{y) for K + 2 < £ < N. 



Since the QNL energy is consistent (see the consistency error analysis in Section [5|), y_F is still an 
equilibrium of £ qnl (y) |25| . Therefore, we will focus on (5 2 £ qnl (y F )u, u) to estimate the stability. 
The second variation of £ qnl (y) evaluated at y = is given by 



{5 2 £ qnl (y F )u,u) = {5 2 £ qnl (y F )u,u) + (5 2 £ qnl (y F )u, u). 



(4.15) 



We first compute the second term of (|4.15|) and get 



(5 2 £ qnl (y F )u,u) 

K 
l=-K 



+ <t>2F K + u 'e-i) 2 + [ u e+i + U 'e+2Y 



]} 



l) 2 + { u 'k+2) 2 +4>2F (%+l + %) 2 + 4 (%+2) 2 } (4.16) 



+ 



2 + (?4+3) 



+ 4>2F 



+ u' K+ i) 2 + 4(^+3) 2 } 



+ ••• 



l \4>"f [{u'k+2, 
N 

£ \ {<& [«)' + (^ + 2 ] + <&F [4 K) 2 + 4 (u' e+1 ) 2 } } ■ 



e=K+3 ' 
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Here we omit the terms whose indices I G {—N + 1, . . . , —{K + 3)} since the QNL energy is 
symmetric. Then we compute the first term, which is given by the following expression: 



(5 2 ^ nt (y F )u,u) = ... 

K 

+ e ^2{ G F [PF( u 'e + u 'e+i) + P2F(u'e-i + u' e + u 'e+i + u 'e+2 )] ' 



+ G' F [p" F {u' t ) 2 + p2 F (u' e + ?4_i) 2 + p' F (u' e+1 ) 2 + ^VK+i + %+2) 21 
+ 2eG" F [p F u' K+1 + p' 2F (u' K+1 + u' K )] 2 + eG' F p" F {u' K+l ) 2 + p 2F (u' K+1 + u' K )' 
+ 2eG" F (p' F + 2p' 2F ) 2 (u' K+2 ) 2 + eG' F (p F + 4p 2F ) (u' K+2 ) 2 

+ 2eG" F [p' F u' K+2 + p' 2F (u' K+2 + u' K+1 )] 2 + eG' F p" F {u' K+2 f + f% F (u' K+2 + u' K+1 
+ 2eG" F (p' F + 2p' 2F f (u' K+3 ) 2 + eG' F {p" F + 4^) (*4+a) 2 

+ e £ [2C; + 2p' 2F f + G' F (p' F + 4/&p)] [(^) 2 + (4 +1 ) 2 ] . 



(4.17) 



Now we use (|4.5p again to rewrite ()4.17j) in the following form 



A? 



2 + K + i) 2 ] 



<J 2 £ v*(y F )u, u) = 6 J] + 2p' 2F ) 2 + G' F {p" F + Ap>> F )\ {{u' t ) 

e=-N+i 

K 

+ ■ ■ ■ ~ e3 E { G "f ■ [(Pf) 2 + 2 °(P2f) 2 + Wf&f] + G'f ■ 2&f) (d {2) u £ 

- e 3 {G'p ■ [{p'p) 2 + W(p' 2F ) 2 + Up F p' 2F ] + G' F ■ 2p' 2 'p} (p<®UK+i 

- e 3 {G F ■ [8(p' 2F ) 2 + 5p' F p' 2F ] + G'p ■ 2p'i F } (d^u k+2 



K+l 



+ e 5 Y,G'p-[8(p' 2 p) 2 + 2p'pp' 2F ] (d® 



K+2 



+ e 5 G' F ■ [4(p' 2F ) 2 + p' F p' 2F ] (D^u K+2 ) - e 7 £ G'p ■ (p' 2F ) 2 (D^ 



'up 
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Combining (S 2 £ qnl (y F )u, u) and (5 2 £ qnl (y f)^, u) together we obtain 



N 



(S 2 £ qnl (y F )u,u)=e £ (i F + A F (Du,) 2 + . . . 



£=-iV+l 
K 



-6 3 E + • [(/>f) 2 + 20(^ 2F ) 2 + W F p f 2F] + G' F ■ 2p>> F } (DM 



1=0 



- e 3 {4>2F + G" F • [{p'f? + nP2F? + n P ' F p' 2F ] + G' F ■ 2p'i F } (d^uk+i 

- 6 3 {4 F + G" F • [8(p' 2F ) 2 + 5p'fP2f] + G' F ■ 2& F } (d^u k+2 



K+l 



+ e 5 Y,G , F-W 2F ) 2 + 2p'FP2F] (p 



(3) 



U( 



e=o 



K+2 



+ e 5 G F • [A(p' 2F f + p' F p'2F] (D^u K+2 ) - e 7 £ G" f • (p' 2F ) 2 [d^u, 



1=0 



Because of the hypotheses ([4.10p and (|4.13p . we have that 

4>2F + G F ■ [{p' F f + lS{p 2F ? + llp' F p' 2F \ + G' F ■ 2p'i F < 0, 

4>2F + G'p ■ [8{p' 2F ) 2 + 5p' F p' 2F ) + G'p ■ 2p 2F < 0. 

Thus, using 



U( 



- 2 r 



and noting that G F • (p' 2F ) 2 > 0, we have 

e> £ G'p ■ [8(p' 2F f + 2p' F p' 2F ] (D^u, 



t=o 



K+2 



+ e 5 G' F ■ [A{p' 2F f + p' F p' 2F ] (d^u k+2 ) - e 7 £ G" f • {p' 2F f (d^u £ 



t=o 



(4.18) 



K+l 



>e"Y,G" F -W 2F f + 2p' F p' 2F ] (D 



)(3) 



1=0 



+ e>G" F ■ [2{p' 2F f + p' F p' 2F ] (d^u k+2 ) > 0. 



So, except in the case K 6 {iV — 2, ... , iV} when there is no continuum region, it follows that y F 
is stable in the QNL model if and only if A F + A F > 0. 

Now we can give a sharp stability estimate for the QNL model from the above estimates and the 
arguments in [9].[T5] . 

Theorem 4.2. Suppose that K < N — 2 and the hypotheses (I4.10P and (I4.13P hold, then the 
uniform deformation y F is stable in the QNL model if and only if A F + A F > 0. 
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Remark 4.2. The role of the assumption (|4.13p in Theorem 14.21 as in Theorem 14.11 is to give 
a necessary condition for u' £ = sm(e£ir) to be the eigenfunction corresponding to the smallest 
eigenvalue of (5 2 £ qnl (y F )u, u) with respect to the norm ||Du|||2. 

Remark 4.3. From Theorem 14.11 and Theorem 14.21 we conclude that the difference between the 
sharp stability conditions of the fully atomistic and QNL models is of order 0(e 2 ). This result is 
the same as for the pair potential case [7J. 

Remark 4.4. We noted in Remark 14.11 that the assumption (I4.13P is necessary for Theorem 14.11 
We now give an explicit example showing that the uniform deformation can be more stable for 
the EAM-QCL model than for the fully atomistic model when (|4.13p fails. We recall that the 
EAM-QCL model is the EAM-QNL model with no atomistic region, that is, 

N 

We consider the case when 



<&F + G" F [p'p + 2p' 2F ) 2 + G' F 2p' 2 ' F > 0, 
which implies that (|4.13p does not hold since it then follows from (I4.10|) that 



,J 2F + G'p 



p' F ) 2 + 20 (p' 2F ) 2 + \2p' F p 2F + G' F 2p 2F 
2F + G" F {p' F + 2p' 2F f + G' F 2p'i F ] + 8G F (2 (p' 2F ) 2 + p' F p' 2F ) 



> 0. 



We define the oscillatory displacement u by 

u e = (-l) e e/(2V2), 



so 



-1)V(V2), pu||^ = l, u'l = (-l)\V2)/e. 
We then calculate from (|4.2p and (|4.4p that 

(5 2 £ a (y F )u,u) = (5 2 £ a (y F )u,u) + (^ a (y F )Ci,u) 



N 1 

6 ^ G' F 2p"/- + + 44b) \\Du\\% + (-ey 2 V)ll^ (2) u 



~ l|2 



-N+l 



Thus, we obtain that 



G' F 2p" F + (0'; + m F ) - m F = 4>" F + G' F 2p" F . 



inf {5 2 £ a (y F )u, u) < f F + G' F 2p" F . 

ueW\{0},||Du|L 2 =l 



On the other hand, we have that 

inf (S 2 £ qd (y F )u, u) = A F +A F = 4 

Therefore, from (I4.19P we have 



(4.19) 



(4.20) 



%f + G"f {p'f + y<2F) 2 + G' F 2pb +( f>' F+ G' F 2p' F . 



inf 

uew\{o},||Du|| 



(5 2 £i d (y F )u, u) > cf) F + G' F 2p" F > 



inf 

uew\{0}, [|Du|l 



(5 2 £ a (y F )u,u) 
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This inequality indicates that the uniform deformation y F can be unstable for the atomistic model, 
but stable for the EAM-QCL model, when the assumption (|4.19p fails. 

We cannot conclude from this argument, though, that the atomistic model is less stable than 
the EAM-QNL model with a nontrivial atomistic region, i.e., K > 0. To see this, we consider an 
oscillatory displacement u£W with support only in the atomistic region (a similar test function is 
used in 0): 



(-lYe 



£ 



2y/2 ' 

0, otherwise. 



(K-1),...,(K-1), 



Then since = (ii£ — ue~i) /e, we have 



£ = K 

2V2 ,c ~ 
0, otherwise. 



■(K-l) 



We substitute the displacement u into (|4.17p and get 

K-3 f 



7 „ , 1 „ 
jPf + 2P2F 



(4.21) 



(S 2 £ " nl (y F )u, u) = e G 'fPf + ^ G" F - [3 (p' 2F ) 2 + 2 (p' F - p' 2F ) 2 \ + G' F 

e=-(K-2) ^ 

= e2(K -2)G' F p" F + 0{e). 
Similarly, we substitute u into (|4.16j) and get 

K-3 

(5 2 £^ l (y F )u, u) = e £ ~<t>" F + 0(e) = e(K - 2)<$ + O(e). (4.22) 

l=-{K-2) 

Therefore, we obtain that 

(5 2 £^ l (y F )u, u) = (6 2 £i nl (y F )u, u) + (5 2 £« nl (y F )u, u) = e(K - 2)($. + 2G' F p" F ) + 0(e). 
Note that 



N 



U 



/||2 



e £ K) 2 = e(if - 1) + - 
e=-N+i 



Thus, we obtain from the above and (|4,20p that 
(5 2 £i nl (y F )u,u) 



u 



/||2 



b F + 2G' F p" F ) + O 



K 



(5 2 £ a (y F )u,u) 



u 



/||2 



+ 



This indicates that when (|4.19p holds and K is sufficiently large, the EAM-QNL model is also less 
stable than the EAM-QCL model. 
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5. Consistency Error and Convergence of The EAM-QNL Model. 

Setting y qnl = y F + u qnl and y a = y F + u a , where both u qnl and u a belong to U, we define the 
quasicontinuum error to be 

e qnl ._ y.a _ qnl _ u a _ u qnl 

To simplify the error analysis, we consider the linearization of the atomistic equilibrium equations 
(|3.3|) and the associated EAM-QNL equilibrium equations (|3.5p about the uniform deformation y F . 
The linearized atomistic equation is 

- {5 2 £ a (y F ) u a , w) = (6T(y F ),w) for all w G U, (5.1) 

and the linearized EAM-QNL equation is 

- (6 2 £ qnl (y F )u qnl ,w) = {5T(y F ),w) for all w e U. (5.2) 

We thus analyze the linearized error equation 

(5 2 £ qnl (y F ) e qnl ,w) = (T qnl ,w) for all weW, (5.3) 

where the linearized consistency error is given by 

<T^, w> := (5 2 £ qnl (y F ) u a , w) - (5 2 £ a (y F ) u a , w) 

= (5 2 £ qnl (y F ) u a , w) - (6 2 £ a (y F ) u a , w) (5.4) 
+ (5 2 £ qnl (y F ) u a , w) - (5 2 £ a (y F ) u a , w). 

Now we will give an estimate of the consistency error T qnl ,w) in the following theorem. We first 
define 

\\ W \\P C {C) := e Z_^ v t' H v lu|(x) : = e 2^^' anci H V II^°°(X) := max tex Nl) for v G 
eec eel 

where C denotes the continuum region {—N + 1, . . . , —{K + 1)} [J{K + 1, . . . , iV} and I denotes 
the interface {-{K + 7), . . . , -K} \J{K, ...,K + 7}. 

Theorem 5.1. The consistency error T qnl ,w), given in (15. 4j) . satisfies the following negative norm 
estimate 



(T qnl , w) < [e 2 [G" F ■ {{p' F ) 2 + l2p' F p' 2F + 20(p' 2F ) 2 ) - 2G' F ■ p'i F + \<% F \] ■ \\D®u a \\ im 

+e 3/2 (Ci + C 2 ) ||£> (2) u a || W) } ||Dw||^ for all weW. 

Proof We focus on the first term of (15. 4p 

<5 2 £^ (y F ) u a , w) - <5 2 £ a (y F ) u<\ w) = • • • + I + I x + I 2 + I 3 , 

where Io is associated with £ = 0, . . . ,K, Ii is associated with £ = K + 1, I 2 is associated with 
£ = + 2 and I3 is associated with £ = K + 3, . . . , N . 

We first compute I3. Note that u a and w are 2./V-periodic, so in the calculation, when the indices 
£ + i > N, i = 1, 2, we can move these terms to the {— N + 1, . . . , — 1} part by using the periodicity 
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as done in (6.9) in [15]. Hence, we can rearrange the terms in I3 to get 



1(5 



J 3 = e Yl G'f- {p'f) 2 {-Du a e ^ + 2Du a e - Du a e+1 ) D Wi (5.5) 

£=K+5 
N 

+ e G "f- {p'fPzf) [ 4 {-Dut-i + 2Du% - Du a e+1 ) + 2 (-Du^_ 2 + 2Du a e - Du a l+2 )} Dw e 

l=K+5 
N 

+ e ^ G' F - (p' 2F ) 2 [3(-Du a i _ 1 + 2Du a e - Du a i+1 ) +2(-Du a e _ 2 + 2Du a e - Du a e+2 ) 
e=K+5 

+ {-Dut_ 3 + 2Du\ - Du a l+3 )] Dw £ 

N 

+ e 2G 'f ■ P2F {-Dv-e-i + 2D ^ - Du1 +l ) Dw t + I31 

£=K+5 



where I31 consists of the interfacial terms, i.e., £ S {K, . . . ,K + 7}, and is given by the following 
expression 



I31 = eG F { {p' F ) 2 [{Du a K+3 - Du a K+A ) w' K+3 + (-Du a K+3 + 2Du a K+A - Du a K+5 ) < +4 ] 

+ p'fp'if [~ {Du a K+3 + Du a K+A ) w' K+2 + {6Du a K+3 - Du a K+2 - 3Du a K+A - 2Du a K+5 ) w' K+3 

+ (l2Du a K+A - Du a K+2 - 3Du a K+3 - 4Du a K+5 - 2Du a K+6 ) w' K+A ] 
+ {p2f) 2 [~ (Du a K+2 + Du a K+3 + Du a K+A + Du a K+5 ) w' K+2 

+ " Du a K+2 - 2Du a K+A - 2Du a K+5 - Du a K+& ) w' K+3 

+ (l3Du a K+A - Du a K+2 - 2Du a K+3 - 3Du a K+5 - 2Du a K+6 - Du a K+7 ) w' K+A ] } 
+ zG'fP2f{ ~ {Du a K+2 + Du a K+3 ) w' K+2 + {2Du a K+3 - Du a K+2 - Du a K+A ) w' K+3 
+ {5Du a K+A - Du a K+3 - 2Du a K+5 ) w' K+A }. 



Since Io is associated with £ = 0, . . . , K where the QNL and the atomistic models coincide with 
each other, we have Iq = 0. Similarly, by direct computation we get the following expression for 
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the sum of Ii and I 2 

Ii + 1 2 = eG" F { (p' F ) 2 [(Du a K+ i - Du a K+2 ) (w' K+l - w' K+2 )] 

+ PfPif [{Du a K+1 - Du a K+2 ) w' K + (2Du a K+1 - 2Du a K+2 + Du% - Du a K+3 ) w' K+1 
+ {6Du a K+2 - Du a K - 2Du a K+1 - Du a K+3 ) w' K+2 - (Du a K+1 + Du a K+2 ) w' K+3 ] 
+ (p' 2F ) 2 [{Du a K + Du a K+1 - Du a K+2 - Du K+3 ) (w' K + w' K+1 ) 
+ {7Du a K+2 - Du a K - Du a K+l - Du a K+3 ) w' K+2 

- (Du a K + Du a K+1 + Du a K+2 + Du a K+3 ) w' K+3 ] } 
+ eG' F p£ F { (3Du a K+2 - Du a K+3 ) w' K+2 - (Du a K+2 + Du a K+3 ) w' K+3 ) 
+ eG F { {p' F ) 2 [{Du k+2 - Du K+3 ) (w' K+2 - w' K+3 )] (5.6) 
+ p'fp'if [{Du a K+2 - Du a K+3 ) w' K+1 + (2Du a K+2 - 2Du a K+3 + Du a K+1 - Du a K+4 ) w' K+2 

+ (6Du a K+3 - Du a K+l - 2Du a K+2 - Du a K+A ) w' K+3 - (Du a K+2 + Du a K+3 ) w' K+4 ] 
+ {P2F) 2 [(Du a K+1 + Du a K+2 - Du a K+3 - Du a K+A ) (w' K+1 + w' K+2 ) 
+ (7Du a K+3 - Du a K+1 - Du a K+2 - Du a K+A ) w' K+3 

~ (Du a K+1 + Du a K+2 + Du a K+3 + Du a K+4 ) w' K+4 \ } 
+ eG' F p'i F { (3Du a K+3 - Du a K+4 ) w' K+3 - {Du a K+3 + Du a K+4 ) w' K+4 }. 

Note that we can rewrite the second term of the second line of I3 as 

2 {-DuU + 2Du a e - Du a £+2 ) = 2 (-Du a e _ 2 + 2DuU - Du a e ) + 4 {-DuU + 2Du a e - Du a £+1 ) 

+ 2 (-Du a e + 2Du a e+1 - Du a e+2 ) . 

Similarly, we can rewrite the third term of the third line of I3 as 

(-Du1_ 3 + 2Du1 - Du1 +3 ) = (-Du1_ 3 + 2Du a t _ 2 - Du^) + 2 (-Du^_ 2 + 2Du1_ x - Duf) 

+ 3 (-£>«?_! + 2Du\ - Du1 +1 ) + 2 (-Dut + 2Du a l+1 - Du a £+2 ) 
+ + 2Du a e+2 - Du a e+3 ) . 

Then we combine Ii, \ 2 and I3 together and rearrange the interfacial terms, i.e., i 6 {K, . . . , K+7}. 
We find that the coefficients of the interfacial terms Ii + 1 2 + 131 are perfectly matched so that they 
are of order e, thus we obtain the following estimate by the Cauchy-Schwarz inequality: 

(5 2 £i nl (y F ) u a , w) - (5 2 £ a (y F ) u a , w) 

< { [G" F ■ ((Pf) 2 + WfP2F + 20(p 2F ) 2 ) - 2G' F ■ p' 2 ' F \ e 2 • \\D^u a \\ em (5.7) 

+C 1 e-||Z)( 2 V||, 2(x) }|| J Dw|| £2 

where X is the interface: {K, . . . , K + 7}, and C\ is a constant independent of e. We note that 

K+7 K+7 

P 2 ul!(z) = « E \ Di2) <\ ± P (2) u a II^CD E e = ^P (2) « a n 2 ^(x)- 

£=X £=K 
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Thus, we obtain 

(6 2 £ (y F ) u a , w) - (5 2 £ a (y F ) u a , w) 

< | e 2 [Gp • ((z^) 2 + 12p' F/! / 2F + 20(p' 2F ) 2 ) - 2G' F • f% F ] • p( 3 )u°|| £f(c) 

+e 3/2 Ci||£> (2) u a ||^ (x) } • pw||^2. 

We can estimate the pair potential consistency error, (5 2 £ qnl (y F ) u a , w) — (5 2 £ a (y F ) u a , w), by 
considering the above estimate for an embedding energy G(4>) = <p/2 to obtain 

(8 2 £ qnl (y F ) u a , w) - (5 2 £ a (y F ) u a , w) 

< {e 2 |0' 2 ' F lll^ {3) u||^(c) + C 2 e\\D^u a \\ e 2 {x) } [|£>w[|^ 

< {e 2 \$ F \\\D®u\\ mc) + C 2 e 3 / 2 \\D^u a \\ er(x) } \\Dw\\^. 
Therefore, we obtain the following optimal order estimate for the consistency error 



w) 



< 



{5 2 £ qnl (y F ) u a , w) - <5 2 £ a (y F ) u a , w) 



+ 



< 



<<5 2 ^ (y F ) u a , w) - (5 2 £ a (y F ) u a , w) 
{e 2 [G F • ((p' F ) 2 + 12/V2F + 20( /9 ' 2F ) 2 ) - 2G' F ■ p'i F + \4 F \] ■ \\D®n a y AC) 
+e 3/2 (Ci + C 2 ) ||D (2) u a || £ c» (x) } pwH^ for all w£W. 

We can now give the convergence result for the linearized EAM-QNL model. 



□ 



Theorem 5.2. Suppose that A F + A F > 0, where A F and A F are defined in (|4.7p and ()4.3j) . and 
that (|4.10p and ()4.13|) holds. Then the linearized atomistic problem (|5. 1 1) as u;eZZ as i/ie linearized 
QNL approximation (|5.2j) /iave unique solutions, and they satisfy the error estimate 

\\Dy a - Dy qra \\p e = \\Du a - Du qnl \\ e 2 

g2 \G" F ■ \{p' F ? + l2p' F p' 2F + 20(p' 2F ) 2 ) - 2G' F ■ p» F + \cf>« F \] • ||£>( 3 V||, ?(C) 



< 



^4 F + A F 



+ 



e 3 / 2 (C 1 + C 2 )\\D^u a \ 



A F + A F 



Proof. The error estimate for the EAM-QNL model follows from the error equation (15. 3ft . the 
stability estimate in Theorem 14.21 an d the consistency estimate in Theorem 15.11 □ 

6. Conclusion. 

We describe a one-dimensional QNL method for the EAM potential following |25j . and we study 
the stability and convergence of a linearization of the next-nearest neighbor EAM-QNL energy. We 
identify conditions for the pair potential, electron density function, and embedding function so that 
the lattice stability of the atomistic and the EAM-QNL models are asymptotically equal. These 
condition are necessary to guarantee that u' e = sin(e£-7r) is the eigenfunction corresponding to the 
smallest eigenvalue of (5 2 £ a (y F )u, u) with respect to the norm ||Z)u||^a. 
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We then give a negative norm estimate for the consistency error and generalize the conclusions 
in [6] to the EAM case. We compare the equilibria of the atomistic and EAM-QNL models and give 
an optimal order 0(e 3//2 ) error estimate for the l\ norm of the strain in terms of the deformation 
in the continuum region. 
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